#################################################################################
# Filename: 3_round2_appendices.R
# Description: Takes the cleaned round 2 and creates 
# materials in the Appendices, including:  
#   (1) covariate balance tests on the round 2 survey sample 
#   (2) regressions for main outcomes
#   (3) marginal Effects plots
#
# NOTE: The one table that should be covered in this script but isn't is the 
# differential attrition analysis displayed in Table A9. Unfortunately, due to the
# data structure in this survey round, this table was created based on back-end
# information contained in the survey platform, and can't be replicated using 
# the publicly available data here. Please contact the authors if this issue is
# of concern. 
#################################################################################

### Load library dependencies 
library(tidyverse)
library(broom)
library(DeclareDesign)
library(estimatr)
library(rio)
library(ggplot2)
library(gridExtra)
library(interflex)
library(stargazer)

# Import cleaned data from 1_round2_loaddata
dat <- import("Round_2_clean.Rdata") 


#################################################################################
# (1) Balance tests (Appendix B)
#################################################################################

pdf("balance_tests_round2.pdf", width=10, height=3.5)
par(mfrow = c(1,3),
    oma=c(1,1,0, 0) ,
    mar=c(4,1,1, 1) )

plot(1, xlim=c(-0.6, 0.6), ylim=c(1, 7), axes=F,  xlab="")

text(y = 7:1, x= 0.6,par("usr")[1], labels =  c("Age", "Education Level", "Urban Resident", "CCP Member", "Income", "Gender", "Hours Watching TV"), srt = 0, pos = 2, xpd = TRUE, cex=1.9)


test1 <-  t.test(dat$age[dat$clip=="Japan"], dat$age[dat$clip=="placebo"])
test2 <-  t.test(dat$education[dat$clip=="Japan"], dat$education[dat$clip=="placebo"])
test3 <-  t.test(dat$urban.residence[dat$clip=="Japan"], dat$urban.residence[dat$clip=="placebo"])
test4 <-  t.test(dat$party[dat$clip=="Japan"], dat$party[dat$clip=="placebo"])
test5 <-  t.test(dat$income[dat$clip=="Japan"], dat$income[dat$clip=="placebo"])
test6 <-  t.test(dat$gender[dat$clip=="Japan"], dat$gender[dat$clip=="placebo"])
test7 <-  t.test(dat$Q14[dat$clip=="Japan"], dat$Q14[dat$clip=="placebo"])


plot(100, pch=16, xlim=c(-0.5, 0.5), ylim=c(1, 7), col="black", xlab="Anti-Japan TV drama", ylab="", yaxt= "n", cex.lab=1.6) 
abline(v=0, col="darkgray", lwd=2)

segments(test1$conf.int[1], 7, test1$conf.int[2], 7, col= 'black', lwd=2)
segments(test2$conf.int[1], 6, test2$conf.int[2], 6, col= 'black', lwd=2)
segments(test3$conf.int[1], 5, test3$conf.int[2], 5, col= 'black', lwd=2)
segments(test4$conf.int[1],4, test4$conf.int[2], 4,  col= 'black', lwd=2)
segments(test5$conf.int[1], 3, test5$conf.int[2], 3, col= 'black', lwd=2)
segments(test6$conf.int[1], 2, test6$conf.int[2], 2, col= 'black', lwd=2)
segments(test7$conf.int[1], 1, test7$conf.int[2], 1, col= 'black', lwd=1)

points(test1$estimate[1]-test1$estimate[2], 7, pch=16, col= 'black', cex=2)
points(test2$estimate[1]-test2$estimate[2], 6, pch=16, col= 'black', cex=2)
points(test3$estimate[1]-test3$estimate[2], 5, pch=16, col= 'black', cex=2)
points(test4$estimate[1]-test4$estimate[2], 4, pch=16, col= 'black', cex=2)
points(test5$estimate[1]-test5$estimate[2], 3, pch=16, col= 'black', cex=2)
points(test6$estimate[1]-test6$estimate[2], 2, pch=16, col= 'black', cex=2)
points(test7$estimate[1]-test7$estimate[2], 1, pch=16, col= 'black', cex=2)


test1 <-  t.test(dat$age[dat$clip=="Anti-American\nSocial Media Video"], dat$age[dat$clip=="placebo"])
test2 <-  t.test(dat$education[dat$clip=="Anti-American\nSocial Media Video"], dat$education[dat$clip=="placebo"])
test3 <-  t.test(dat$urban.residence[dat$clip=="Anti-American\nSocial Media Video"], dat$urban.residence[dat$clip=="placebo"])
test4 <-  t.test(dat$party[dat$clip=="Anti-American\nSocial Media Video"], dat$party[dat$clip=="placebo"])
test5 <-  t.test(dat$income[dat$clip=="Anti-American\nSocial Media Video"], dat$income[dat$clip=="placebo"])
test6 <-  t.test(dat$gender[dat$clip=="Anti-American\nSocial Media Video"], dat$gender[dat$clip=="placebo"])
test7 <-  t.test(dat$Q14[dat$clip=="Anti-American\nSocial Media Video"], dat$Q14[dat$clip=="placebo"])


plot(100, pch=16, xlim=c(-0.5, 0.5), ylim=c(1, 7), col="black", xlab="Anti-US social media clip", ylab="", yaxt= "n", cex.lab=1.6) 
abline(v=0, col="darkgray", lwd=2)

segments(test1$conf.int[1], 7, test1$conf.int[2], 7, col= 'black', lwd=2)
segments(test2$conf.int[1], 6, test2$conf.int[2], 6, col= 'black', lwd=2)
segments(test3$conf.int[1], 5, test3$conf.int[2], 5, col= 'black', lwd=2)
segments(test4$conf.int[1],4, test4$conf.int[2], 4,  col= 'black', lwd=2)
segments(test5$conf.int[1], 3, test5$conf.int[2], 3, col= 'black', lwd=2)
segments(test6$conf.int[1], 2, test6$conf.int[2], 2, col= 'black', lwd=2)
segments(test7$conf.int[1], 1, test7$conf.int[2], 1, col= 'black', lwd=2)

points(test1$estimate[1]-test1$estimate[2], 7, pch=16, col= 'black', cex=2)
points(test2$estimate[1]-test2$estimate[2], 6, pch=16, col= 'black', cex=2)
points(test3$estimate[1]-test3$estimate[2], 5, pch=16, col= 'black', cex=2)
points(test4$estimate[1]-test4$estimate[2], 4, pch=16, col= 'black', cex=2)
points(test5$estimate[1]-test5$estimate[2], 3, pch=16, col= 'black', cex=2)
points(test6$estimate[1]-test6$estimate[2], 2, pch=16, col= 'black', cex=2)
points(test7$estimate[1]-test7$estimate[2], 1, pch=16, col= 'black', cex=2)

dev.off()


#################################################################################
# (2) Regressions for main outcomes
#################################################################################
# Table A3 - Round Two: OLS Results for Main Outcomes
test1 <- lm(Q19_1~clip, data=dat )
test2 <- lm(Q19_1~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )
test3 <- lm(scale(antiforeign_usa_pca)~clip, data=dat )
test4 <- lm(scale(antiforeign_usa_pca)~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )
test5 <- lm(scale(antiforeign_japan_pca)~clip, data=dat )
test6 <- lm(scale(antiforeign_japan_pca)~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )

summary(test1)
stargazer(test1, test2, test3, test4, test5, test6,
          covariate.labels = c("Anti-Japan television drama", "Anti-US social media clip", "Age", "Education", "Urban resident", "Party", "Income", "Gender", "Hours watching TV"),
          dep.var.labels = c("Anger", "Anti-American sentiment", "Anti-Japanese sentiment"),
          omit.stat=c("f", "ser"))


#Table A4 - Round Two and Three: OLS Results for Main Outcomes
test1 <- lm( protest.gov~clip, data=dat )
test2 <- lm(protest.gov ~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )
test3 <- lm( protest.gov_f~clip, data=dat )
test4 <- lm(protest.gov_f ~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )
test5 <- lm(responsive.gov ~clip, data=dat )
test6 <- lm( responsive.gov~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )

summary(test1)
stargazer(test1, test2, test3, test4, test5, test6,
          covariate.labels = c("Anti-Japan television drama", "Anti-US social media clip", "Age", "Education", "Urban resident", "Party", "Income", "Gender", "Hours watching TV"),
          dep.var.labels = c("Would protest government (Wave 2)", "Would protest government (Wave 3)", "Government is sensitive"),
          omit.stat=c("f", "ser"))


# Table A5 - Round Three: OLS Results for Main Outcomes
test1 <- lm(Q19_1_f~clip, data=dat )
test2 <- lm(Q19_1_f~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )
test3 <- lm(scale(antiforeign_usa_pca_f)~clip, data=dat )
test4 <- lm(scale(antiforeign_usa_pca_f)~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )
test5 <- lm(scale(antiforeign_japan_pca_f)~clip, data=dat )
test6 <- lm(scale(antiforeign_japan_pca_f)~clip+age+education+urban.residence+party+income+gender+Q14, data=dat )

summary(test1)
stargazer(test1, test2, test3, test4, test5, test6,
          covariate.labels = c("Anti-Japan television drama", "Anti-US social media clip", "Age", "Education", "Urban resident", "Party", "Income", "Gender", "Hours watching TV"),
          dep.var.labels = c("Anger", "Anti-American sentiment", "Anti-Japanese sentiment"),
          omit.stat=c("f", "ser"))




###Multiple Comparisons Testing####
lm1 <- lm(Q19_1~clip, data=dat )
lm2 <- lm(antiforeign_usa_pca~clip, data=dat )
lm3 <- lm(antiforeign_japan_pca~clip, data=dat )
lm4 <- lm(militant_pca~clip, data=dat )
lm5 <- lm(patriot_pca~clip, data=dat )
lm6 <- lm(responsive.gov~clip, data=dat )
lm7 <- lm(protest.gov~clip, data=dat )

p <- c(summary(lm1)$coefficients[2,4],
       summary(lm1)$coefficients[3,4],
       summary(lm2)$coefficients[2,4],
       summary(lm2)$coefficients[3,4],
       summary(lm3)$coefficients[2,4],
       summary(lm3)$coefficients[3,4],
       summary(lm4)$coefficients[2,4],
       summary(lm4)$coefficients[3,4],
       summary(lm5)$coefficients[2,4],
       summary(lm5)$coefficients[3,4],
       summary(lm6)$coefficients[2,4],
       summary(lm6)$coefficients[3,4],
       summary(lm7)$coefficients[2,4],
       summary(lm7)$coefficients[3,4])

p.adjust(p, method = "BH", n = length(p))


lm1 <- lm(Q19_1_f~clip, data=dat )
lm2 <- lm(antiforeign_usa_pca_f~clip, data=dat )
lm3 <- lm(antiforeign_japan_pca_f~clip, data=dat )
lm4 <- lm(militant_pca_f~clip, data=dat )
lm5 <- lm(patriot_pca_f~clip, data=dat )
lm6 <- lm(responsive.gov_f~clip, data=dat )
lm7 <- lm(protest.gov_f~clip, data=dat )

p <- c(summary(lm1)$coefficients[2,4],
       summary(lm1)$coefficients[3,4],
       summary(lm2)$coefficients[2,4],
       summary(lm2)$coefficients[3,4],
       summary(lm3)$coefficients[2,4],
       summary(lm3)$coefficients[3,4],
       summary(lm4)$coefficients[2,4],
       summary(lm4)$coefficients[3,4],
       summary(lm5)$coefficients[2,4],
       summary(lm5)$coefficients[3,4],
       summary(lm6)$coefficients[2,4],
       summary(lm6)$coefficients[3,4],
       summary(lm7)$coefficients[2,4],
       summary(lm7)$coefficients[3,4])

p.adjust(p, method = "BH", n = length(p))


#################################################################################
# (3) Marginal effects plots
#################################################################################
dat$age_graphics <- ifelse(dat$age==1, 18, dat$age*5 + 15 )

pdf("r2_age_anger_jp_ME.pdf")
inter.binning(Y = "Q19_1", D = "japan", X = "age_graphics", data = dat, nbins=5, 
              Xlabel="Age", Dlabel = "Anti-japanese \n television drama", Ylabel="anger",
              show.grid=F,theme.bw=T, vartype = "robust", main = NULL, na.rm=T)$graph
dev.off()

pdf("r2_age_anger_hk_ME.pdf")
inter.binning(Y = "Q19_1", D = "hk", X = "age_graphics", data = dat, nbins=5, 
              Xlabel="Age", Dlabel = "American interference\n in Hong Kong video", Ylabel="anger",
              show.grid=F,theme.bw=T, vartype = "robust", main = NULL, na.rm=T)$graph
dev.off()

pdf("r2_age_antif_jp_ME.pdf")  
inter.binning(Y = "antiforeign_japan_pca", D = "japan", X = "age_graphics", data = dat, nbins=5, 
              Xlabel="Age", Dlabel = "Anti-japanese television\n drama", Ylabel="Anti-Japanese sentiment",
              show.grid=F,theme.bw=T, vartype = "robust", main = NULL, na.rm=T)$graph
dev.off()

pdf("r2_age_antif_hk_ME.pdf")   
inter.binning(Y = "antiforeign_hk_pca", D = "hk", X = "age_graphics", data = dat, nbins=5, 
              Xlabel="Age", Dlabel = "American interference in Hong\n Kong video", Ylabel="Anti-American sentiment",
              show.grid=F,theme.bw=T, vartype = "robust", main = NULL, na.rm=T)$graph
dev.off()





